Weak-coupling quantum Monte Carlo calculations on the Keldysh contour: theory and application 
to the current- voltage characteristics of the Anderson model 
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We present optimized implementations of the weak-coupling continuous-time Monte Carlo method defined 
for nonequilibrium problems on the Keldysh contour. We describe and compare two methods of preparing the 
system before beginning the real-time calculation: the "interaction quench" and the "voltage quench", which are 
found to be suitable for large and small voltage biasses, respectively. We also discuss technical optimizations 
which increase the efficiency of the real-time measurements. The methods allow the accurate simulation of 
transport through quantum dots over wider interaction ranges and longer times than have heretofore been pos- 
sible. The current-voltage characteristics of the particle-hole symmetric Anderson impurity model is presented 
for interactions U up to 10 times the intrinsic level width T. We compare the Monte Carlo results to fourth order 
perturbation theory, finding that perturbation theory begins to fail at U/Y > 4. Within the parameter range 
studied we find no evidence for a splitting of the Kondo resonance due to the applied voltage. The interplay of 
voltage and temperature and the Coulomb blockade conductance regime are studied. 
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I. INTRODUCTION 

The development of robust methods for the computation of 
nonequilibrium properties of quantum many-particle systems 
is a crucial issue in present-day condensed matter physics, 
with impact on topics ranging from nonequilibrium transport 
in nanostructures 1 to pump-probe spectroscopy of bulk con- 
densed matter systems 2 ' 3 and the wide range of new spectro- 
scopies possible in cold atom systems. 4 An important step 
forward occurred with the development of continuous-time 
quantum Monte Carlo (CTQMC) methods for impurity mod- 
els. These algorithms were first introduced as imaginary-time 
methods for obtaining equilibrium properties 5 ' 6 ' 7,8 and soon 
afterwards were extended to real-time dynamics and nonequi- 
librium problems. 9 ' 10 ' 11 ' 12 

The continuous-time methods are in essence stochastic 
samplings of diagrammatic expansions of the time evolution 
operator. The mean perturbation order required in the cal- 
culation increases with the time (or inverse temperature) to 
be studied and the calculations are limited by the perturba- 
tion order which can be achieved with given computational re- 
sources. In the equilibrium case one considers the imaginary- 
time evolution operator exp[— tH] which is real and positive 
definite, so the computational task is to estimate a sum of 
real (decaying) exponentials and the only sign problem which 
arises is the fermion sign problem occurring in models com- 
plicated enough to sustain fermion loops. For these reasons 
the CTQMC methods have proven to be very powerful in the 
equilibrium context. 13 In the real time case, on the other hand, 
one must consider the intrinsically complex time evolution 
operator exp[— itH], and convergence comes from the can- 
cellation of oscillations. The theoretical task is therefore to 
estimate the sum of terms with oscillating signs or rotating 
phases and a severe "dynamical" sign problem occurs even 
in the absence of fermion loops. The average sign decreases 
exponentially with perturbation order, which limits the acces- 



sible range of interaction strengths and simulation times. 

Because of these limitations, important questions such as 
the nonequilibrium Kondo effect could so far not be ade- 
quately addressed. The equilibrium Kondo effect in quan- 
tum dots, 14,15 which involves the formation of a scattering 
resonance (density of states peak) at the Fermi level, was 
experimentally confirmed in the zero bias limit. 1 While at 
very low voltage biasses, the pinning of the Kondo reso- 
nance to the Fermi level leads to an unrenormalized conduc- 
tance for symmetric dots, it is well known that a high volt- 
age bias destroys the Kondo effect. The crossover from the 
low voltage universal regime ("linear response regime") to the 
higher bias Coulomb blockade regime is presently not under- 
stood. It has been proposed on the basis of the noncrossing 
approximation, 16 real time diagrammatic methods, 17 and per- 
turbative calculations 18 that the peak in the density of states 
splits into two in a certain parameter regime. Our previous 
investigation of the non-equilibrium Anderson model 11 pro- 
duced no sign of this phenomenon, but the accuracy of the 
simulations in the low-bias region was not sufficient to set- 
tle the issue. Methodological improvements allowing a more 
accurate numerical study of the small-to-intermediate voltage 
regime are therefore needed. 

The existing continuous-time Quantum Monte Carlo ap- 
proaches for nonequilibrium systems are more-or-less direct 
extensions of the imaginary time algorithms previously devel- 
oped. It appears worthwhile to attempt to optimize them, even 
though the dynamical sign problem inherent in these methods 
unavoidably limits what can be achieved. In this paper we 
present an efficient implementation of the weak-coupling dia- 
grammatic Monte Carlo method for non-equilibrium systems, 
describing ways to reorganize the expansion and to improve 
the measurement formulae in order to increase the accuracy 
of the numerical data for a given set of parameters. 

The method introduced previously 11 corresponds to the 
simulation of a system prepared in the nonequilibrium but 
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non-interacting state, with the interaction turned on at time 
t = 0. We refer to this simulation method as an "interaction 
quench". Since the real-time methods compute the time evolu- 
tion of the system after the quench, an important consideration 
is the time needed for the system to evolve to the interacting 
steady state. Optimized preparation of the initial ensemble has 
the potential to reduce this relaxation time, therefore leading 
to simulations requiring a smaller total time interval for the 
measurement of a given property. Motivated by this idea we 
extend the formalism from two real-time branches to an "re- 
shaped" contour which includes an imaginary time branch. 
Evolution along the imaginary time branch may be thought of 
as preparing the system in a correlated equilibrium state, after 
which the voltage is turned on at time t = 0. We refer to this 
simulation method as a "voltage quench". 

One purpose of this paper is to compare interaction and 
voltage quenches. We will show that at temperature T = in- 
teraction quenches are suitable for voltage biasses larger than 
the Kondo temperature (i. e. for voltage biasses large enough 
to suppress the Kondo resonance in the many-body density 
of states). The times which can be reached in the Monte 
Carlo simulation are long enough to observe convergence into 
a steady state even at large interaction strengths. On the other 
hand, if the voltage bias is small and the temperature is finite, 
the voltage quench is a suitable alternative, because it allows 
the important ground state correlations to be built up via the 
computationally less problematic imaginary time evolution. 

We show that our optimized implementation allows the 
computation of accurate currents over a wide voltage range, 
even for interaction strengths which are clearly outside the 
reach of low-order perturbation theory. We use the numeri- 
cal results to test predictions based on fourth order perturba- 
tion theory in the self-energy. 18 We determine the largest in- 
teraction strength for which the perturbation theory provides 
accurate results over the entire voltage range, and for larger in- 
teractions, the voltage window where deviations appear. The 
predicted splitting of the Kondo resonance 161718 is not evi- 
dent in the numerical data. 

The rest of this paper is organized as follows. In section II 
we introduce the model to be solved and present the methods 
used to solve it, in particular defining the voltage and interac- 
tion quenches. Sections III and IV present results for the in- 
teraction and voltage quenches respectively. Section V gives 
results for the current-voltage characteristics of the model and 
section VI is a summary and conclusion. 



the terms 

a=L,R p,<y 

H mlx = E E(W^ + ^ C -), ( 2 ) 

a=L,R p.o 

H dot = £d^2n d . a , (3) 

a 

H L r = U(n d ^n dtl -(n dA + n dtl )/2). (4) 

In the following we will consider two sources of time depen- 
dence in Hqj. In the interaction quench we take U = for 
times t < with an instantaneous step to a non-zero U at 
t = 0; in the voltage quench we take \xl = \ir for time t < 
with an instantaneous step to a nonzero /i£ — at t = 0. We 
assume that the lead electrons equilibrate instantly to the new 
chemical potential so that the equal time correlators of lead 
operators are (a^ a a^ a ,) = 5 a ^5 p ^5 a y frA^a - Ma), 
with ./t(-c) = {e x ^ T + the Fermi distribution function 
for temperature T and fj, a the value of the chemical potential 
for lead a at the appropriate time. 

In this paper we will consider only symmetric voltage bi- 
asses (hl = — lift = V/2) and half-filled dots (e d — 0). The 
consequences of relaxing these assumptions will be briefly 
mentioned in the conclusions. The energy scales of the model 
are set by the level broadenings 

r° =7rXX<f (5) 
p 

associated with the leads a. The total level broadening 

r = r L + r R (6) 

is used as the energy unit throughout the paper. 

We consider fiat bands centered at zero, with a high en- 
ergy cutoff u c . As we shall see in Section IV a sharp high 
frequency cutoff leads to oscillations in the time evolution of 
the current after a voltage quench. A sufficiently smooth band 
cutoff damps the oscillations but does not affect the steady 
state current. We adopt a Fermi-function like smoothing with 
"smoothing parameter" v, 

yL,R 

Y L ' R (uj) = - (7) 



II. MODEL AND METHODS 
A. Model 

We consider the one-orbital Anderson impurity model, 
which describes a single spin-degenerate (<r) level with a Hub- 
bard interaction U (the "dot") coupled by hybridization V to 
two reservoirs ("leads") labeled by a = L, R. The Hamilto- 
nian Hqj = H® ot + Hjj + i?bath + -^mix of this model contains 



B. Real-time Monte Carlo method: weak-coupling approach 

We use the weak-coupling formulation of the real-time di- 
agrammatic Monte Carlo approach as described in Ref. 11. 
This is a real-time implementation of the continuous-time 
auxiliary field algorithm 8 which is based on the combination 
of a weak-coupling expansion and an auxiliary field decom- 
position. Here, we will briefly summarize the main aspects 
of this method and then discuss some relevant issues concern- 
ing its efficient implementation. In order to enable simula- 
tions starting from an interacting initial state, we formulate 
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FIG. 1 : Illustration of the Keldysh contour for the interaction quench 
(top panel) and voltage quench (bottom panel). In an interaction 
quench starting from U = 0, the imaginary time branch of the con- 
tour is shifted to t — —00 and need not be explicitly considered in 
the Monte Carlo simulation. The red arrows represent auxiliary Ising 
spin variables. The top panel shows a Monte Carlo configuration cor- 
responding to perturbation order n+ = 2, n_ = 2, and the bottom 
panel a configuration corresponding to n+ = 3, n_ = 2, np = 2. 



with Kp and K t some arbitrary (non-zero) constants. The in- 
teraction and the chemical potentials need not be the same 
on the imaginary time branch as they are on the real-time 
branches. In the formalism as written the interaction and 
chemical potentials are taken to be time independent on the 
real-time branches, but it is straightforward to generalize the 
method to time-dependent U and fi. The notation indi- 
cates that on the real-time portion of the contour the two leads 
have different chemical potentials fj, a = fj,o ± Sfi, whereas 
iJ^ath means tnat on tne imaginary time portion of the contour 
the two leads have the same chemical potential /xn. Hence- 
forth we choose energies such that jUrj = and consider a 
symmetrically applied bias voltage V (<5/i = V/2). 

The time evolution along the real-time and imaginary-time 
contours is expanded in powers of Hjj — K t /t and Hjj — 
Kp I (3, respectively. Each interaction vertex is then decoupled 
using Ising spin variables according to the formula 5 (x = t or 



Hu-Kx/x = e 

s=-l,l 

cosh( 7x ) = 1 + {xU)/(2K 



7s(™<J,T-™d,l) ^ (9) 

(10) 



the method on the L-shaped contour which runs from to t 
and back to along the Keldysh real time axis, then to — i/3 
along the imaginary time axis. 

The weak coupling algorithm may be taken to start from the 
following expression for the partition function Z = Tre~° H : 

Z = e K eTr[e- 0(H ^ +H ^ +H ^ +H o- K e /l3) 

xe *(H'Z+ H L+ H u+H m ,-K t /t) 
xe -it(H^ 1 +H°, t +Hu+H mix -K t /t)l ^ (g) 



The resulting collection of Ising spin variables on 
the contour represents the Monte Carlo configuration 
{(<i, si), (t 2 ,s 2 ), ■ ■ ■ (t n ,s n )}, with U denoting the position 
of spin i on the L-shaped contour (see illustration in Fig. 1). 
There are n + spins on the forward branch, n_ spins on the 
backward branch and np spins on the imaginary-time branch 
of the contour (n = n + + n_ + np). The weight of such 
a configuration is obtained by tracing over the dot and lead 
degrees of freedom and can be expressed in terms of two 
determinants ofnxn matrices iV^ 1 : 8 



ti>({(fi,si),(t2,*2),... (*»,*»)}) = H"-) (t n+ )(K t dt/2t) n - +n +(K p dr/2f3)^l[detN- 1 , (11) 

(T 

N- 1 = e s ° ~ (iG ,„)(e s ° - I). (12) 
I 



Here (G"o,o-)ij = Go.<r(U,tj) is the ij element of the n x n 
matrix of non-interacting Green functions 

G Q!(T (t,t') = -i(Tcd a (t)dt(t')) (13) 

computed using the possibly time-dependent chemical poten- 
tials and evaluated at the time arguments defined by the Ising 
spins. The quantity e s ° = diag(e 7lSlcr , . . . , e 7 " SntT ) is a di- 
agonal matrix depending on the spin variables (with 7; = 74 
for spins located on the real-time branches and 7, = 7 ^ for 



spins on the imaginary time branch). Tc is the contour order- 
ing operator, which exchanges the product A(t)B(t') of two 
operators if t is earlier on the contour than t' (a minus sign is 
added if the exchange involves an odd number of Fermi oper- 
ators). 

A Monte Carlo sampling of all possible spin configura- 
tions is then implemented based on the absolute value of the 
weights (11). The contribution of a specific configuration 
c = {(*i,si), (£2,^2), • ■ ■ (t n , s n )} to the current is given by 11 
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A c a (t,t') = A 0>rT (t,t')+i ]T G 0!<7 (t,U)[(e s " - ^N^ijAo^t^t'), (14) 

I 



with the first term on the right hand giving the contribution to 
the non-interacting current and the second term a correction 
due to the interactions. In Eq. (14) 

Ao, a (t,t') = {T c ~a L a \t')d a {t)) Q (15) 

denotes a dot-lead correlation function of the noninteracting 
model for the composite left lead operator = Vp a pa- 
The current expectation value is 

I(t) = -2Im^[<^(M)0 c )/(0 c )], (16) 



Go(M') 
A (t,t') 



In the voltage quench, on the other hand, the interaction is 
non-vanishing on the imaginary time portion of the contour 
(Fig. 1), while the chemical potential difference jumps instan- 
taneously from zero (on the imaginary branch) to V (on the 
real branches). Because of the time dependence of the chemi- 
cal potentials, the noninteracting Green functions are not time 
translation invariant and we cannot express Go,<j and the dot- 
lead correlator Ao, a in the form of a Fourier transform. In- 
stead, those functions must be computed numerically from 
their equations of motion, as explained in the appendix. 

C. Optimization of the Monte Carlo sampling 

The sign (phase) problem in the weak-coupling CTQMC 
method grows exponentially with the average perturbation or- 
der on the real-time branches, which in turn is proportional to 
the simulation time, while operators on the imaginary time 
branch do not add significantly to the sign problem. To 
reach long times or strong interactions, it is therefore impor- 
tant to reduce the average perturbation order on the real-time 
branches as much as possible. An essential point to note in 



where (.) denotes the Monte Carlo average and <fi c the phase 
of the weight of the configuration c. 

In an interaction quench, the imaginary-time evolution is 
not explicitly considered in the Monte Carlo simulation and 
temperature appears only as a parameter in the noninteracting 
Green functions (see Fig. 1). Moreover, the latter depend only 
on time differences, and thus can be easily expressed in terms 
of their Fourier transform. Assuming a large band cutoff and 
neglecting the real part of the lead self-energy we find 1119 



(17) 



(18) 



this context is that in the particle-hole symmetric case, the pa- 
rameters K x of the algorithm can be chosen such that only 
even perturbation orders appear in the expansion. In fact, for 

Kx = -xU/4 (19) 
the spin degree of freedom effectively disappears (e 7Scr = — 1) 
and the algorithm becomes the real-time version of Rubtsov's 
weak-coupling method 6 for the particle-hole symmetric inter- 
action term Hjj — K x /x = U{n d ,] — \){n d ,i — \)- (For 
a detailed discussion of the equivalence between the Rubtsov 
and CTAUX methods for the Anderson impurity model, con- 
sult Ref. 21). The odd perturbation orders are continuously 
suppressed as K x approaches — xU /4. For K x = —xU/4 + 5 
and sufficiently small <5, the average perturbation order can 
be reduced by about half compared to the \K\ = 0.1 used in 
the simulations presented in Ref. 11. This in turn allows us 
to reach times and interaction strengths which are a factor of 
two larger. We note in passing that the suppression of odd 
perturbation orders was also essential in the nonequilibrium 
dynamical mean field calculations of Ref. 20. 

We next discuss some tricks to improve the efficiency of the 
current measurement. First, we rewrite Eq. (14) as 



v fdw _^ r"M(/(a,- Ata )-e c (M')) 

Z f n J 2n e ( W -e d - C//2)2 + ' 

■ f &± -Mt-t') r L (u)r R (u){f(uj - hl) - /(cj - n R )) 

J 2tt (w - e d - L//2) 2 + IV) 2 

(dw _ Mt _ tl) r £ H(q, -e d - U/2)(f{u - » L ) - O c (t, t')) 
J 2ir (w -e d - (7/2)2 + r 2 



I 

A c a (t,t') =4>,«r(t,t / ) + / ds i J ds 2 G . a (t, S c(si,ti)i(e S ^ I)N a ] tJ 5 c ( s 2,tj))AoAs2,t'), (20) 
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where the variables si and s 2 run over the entire contour and 
the contour delta function is defined by J ds6c(t, s)f(s) = 
f(t). It is therefore sufficient to accumulate the quantity 

n 

X a { Sll s 2 ) = (* £ fc(»i,*i)[(e s ' -lWijScfatj)). 

(21) 

Furthermore, it follows from Eq. (11) that the weight of a 
Monte Carlo configuration changes sign if the last spin (cor- 
responding to the largest time argument) is shifted from the 
forward contour to the backward contour or vice versa. Since 
the absolute value of the weight does not change, these two 
configurations will be generated with equal probability. As a 
result, all the terms in Eq. (20) which do not involve the last 
operator on the contour will cancel. It is therefore more effi- 
cient and accurate to accumulate 

n 

X rT (s 1 ,s 2 ) = -S({ti})) ^2 x ( s i^> s 2,j) 



id({ti}) X] [^I'last; s 2, + x(si, I; S2,last)] 



l not last 



(22) 



withx(s 1> i;s 2 ,j) = 6 c (si, ii)[(e r " - l)N a \ ltj 5 c {s 2 , t) and 
S({U}) = 1 if maxi Re(tj) > and otherwise. 

Also, by comparing the contributions to the current of 
the original configuration and the one with the last operator 
shifted from the upper to the lower contour (or vice versa), 
one finds that they almost (but not completely) cancel. The 
errorbars on the current can thus be substantially reduced by 
appropriate symmetrizations of X(si, s 2 ). 



III. RESULTS: INTERACTION QUENCH 

A. Convergence to the long time limit: large bias voltage 

Calculations based on an interaction quench from [7 = 
are particularly simple, because there are no interaction ver- 
tices (or spins) on the imaginary time branch, and only the 
real-time branches of the contour need to be considered in the 
simulation. Temperature enters only as a parameter in the lead 
correlators, making it possible to treat arbitrary temperatures, 
including T = 0. At time t — 0, the system is non-interacting 
but subject to an applied bias V, so a current Iq(V) appropri- 
ate to the non-interacting model is flowing through the dot. 
At t = 0+ the interaction is turned on and the system relaxes 
into the steady-state configuration appropriate to the interact- 
ing model. Figure 2 shows the time dependence of the current 
calculated for the large bias voltage V/Y = 4 and several in- 
teraction strengths. We see that the transient behavior is such 
that the current initially decreases sharply, overshoots and 
eventually relaxes more slowly back up into the new steady 
state. The interaction-dependence of the steady-state current 
is a consequence of the Coulomb blockade physics, apparent 
even at the large voltages studied here. 




FIG. 2: Time evolution of the current for V/Y = 4 and different 
interaction strengths (T = 0). In the initial state, the current is given 
by the steady state current through the non-interacting dot. At time 
t = 0, the interaction is turned on. After a time of a few inverse 
r, the current saturates at the value corresponding to the steady state 
current in the interacting dot. 




FIG. 3: Time evolution of the current for different voltage biasses and 
interaction strength U/Y = 6 (T = 0). In the initial state, the current 
is given by the steady state current through the non-interacting dot. 
At time t = 0, the interaction is turned on. 



For intermediate and large voltage bias (V/Y > 2) and not 
too large interaction (U/Y < 8) the time required for conver- 
gence to the steady state is tY w 2, essentially independent 
of interaction strength. Given the scaling of the perturbation 
order (and hence the sign problem) with U and t, interactions 
up to U/Y < 10 are accessible with the current implemen- 
tation. A comparison of Fig. 2 to Fig. 13 of Ref. 11 shows 
that the technical improvements introduced in this paper have 
substantially extended the range of applicability of the weak- 
coupling Monte Carlo method (about a factor 2-3 in U or t) 
and allow us to obtain accurate results in the intermediate-to- 
strong correlation regime. 

In Fig. 3 we plot the time evolution of the current for 
fixed U/Y = 6 and several voltage biasses. For voltages 
V/Y > 2, even though the transient behavior is clearly 
voltage-dependent, the current settles into the new steady state 



after a time tT ps 2. However, as the voltage is decreased be- 
low V/T « 2 the transient time increases. At V = T the 
long time limit is attained only for tT > 3 and as V is further 
decreased the approach to the asymptotic behavior becomes 
even slower. 



B. Convergence to the long time limit: small bias voltage 

To better analyse the approach to steady state at small volt- 
ages we present in the upper panel of Fig. 4 the time depen- 
dence of the current for several smaller voltages and two inter- 
action strengths. For better comparison, we plot here the ratio 
I/Iq of the interacting current I to the noninteracting current 
Iq. One sees that as V is decreased or U is increased the evo- 
lution of the current from the post-quench minimum to the 
long-time steady state value takes an increasingly long time. 
Since the longest accessible time is tT « 6 for U/T = 4 and 
tT w 4 for U/T = 6, the accurate measurement of I becomes 
impossible in the small voltage regime. However the short- 
time transient behavior is accessible at all voltages. While the 
ratio (J/Jo)(i) is clearly voltage dependent at higher biasses, 
the data seem to converge as V is reduced to a non-trivial 
curve with a pronounced minimum near an only weakly U - 
dependent time w 1. 

We believe that the increasingly slow convergence as V — > 
is a signature of the Kondo effect, which is characterized 
by an energy scale which becomes exponentially small as 
U increases. After the interaction quench, the Kondo reso- 
nance has to be built up as time progresses, and in the limit 
V — > 0, T — > this requires an increasingly large number of 
interaction vertices and hence an increasingly long simulation 
time. On physical grounds one expects that the time needed 
to evolve into steady state is proportional to the inverse of the 
associated energy scale. 

Empirically, we find that the slow relaxation becomes an 
issue in the linear response regime, where the non-interacting 
and interacting currents are very similar. For V/T > 0.5, 
where the interacting current is substantially smaller than Iq, a 
useful estimate of I seems possible, even though in the voltage 
window up to V/T « 2 a small drift in the current may remain 
up to the longest accessible times. This drift makes it difficult 
to define reliable error bars on I, but it appears unlikely that 
the steady state value will differ from I(t max ) by more than 
the largest deviation in the window [i max /2, t max ], which we 
use as error estimate. For V/T = 0.25, an accurate estimate 
is not possible from the interaction quench procedure, but the 
current is very close to the linear response value, so the un- 
certainty is in fact not that important. This is illustrated in 
the bottom panel of Fig. 4, which compares the Monte Carlo 
data to the noninteracting current and results from fourth order 
perturbation theory. 18 The plot also indicates as hashed region 
the voltage range V/T < 0.4 where accurate measurements 
of the long-time limit become prohibitively difficult at T = 0. 
We will see below that this roughly corresponds to voltages 
smaller than the Kondo temperature Tk ■ 
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FIG. 4: Interaction quench in the small-voltage regime (T = 0): The 
top panel shows the ratio of interacting to noninteracting current for 
U/T = 4 and U/T — 6 and indicated voltage biasses. For V/T < 
0.5 the time needed to reach the steady state grows much beyond the 
largest time accessible in the Monte Carlo simulation. Bottom panel: 
comparison of (7-quench estimates (symbols, red and blue online) 
for the steady-state current to the noninteracting current (thick black 
line) and fourth order perturbation theory 18 (light lines, red and blue 
online). For V/T > 0.5 the Monte Carlo data show the value of 
/(tmax) with errorbar max t6 [ iim / 2l t„ I ] \I(t) - I(t m ^)\ (t mlx T = 6 
for U/T = 4 and t raax r = 4 for U/T = 6). For V/T = 0.25, we 
use I « (J(imax) + ^o)/2 with errorbar of size (To — I{t ma ))/2. 



C. Temperature Dependence 

It is also of interest to examine the temperature dependence 
of the current. The interplay between voltage and tempera- 
ture as the Kondo regime is approached presents an interest- 
ing problem. One expects that as the temperature is increased, 
the Kondo effect gets washed out and the simulations would 
therefore more readily converge even at small bias voltages. 
The temperature dependence of the current calculated from 
the interaction quench for U/T = 6 and several values of the 
voltage bias is plotted in Fig. 5. In the linear response regime 
(V/T = 0.125, 0.25) the ratio of the interacting current I(T) 
to the noninteracting current Iq(T = 0) exhibits a strong tem- 
perature dependence, even at T/V <C 1. The temperature 
dependence arises because lowering the temperature strength- 
ens the Kondo resonance and leads to an increase in the in- 
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T/r 

FIG. 5: Temperature dependence of the ratio of interacting current 
at temperature T to noninteracting current at temperature zero for 
indicated values of the voltage bias. The interaction strength is 
U/T = 6. The symbols show Monte Carlo results, the black line 
the analytical curve for V — * extracted from Ref. 25 and plotted 
forTx/r = 0.24. 



teracting current. The temperature dependence for small volt- 
age bias (V/T = 0.125) approaches the analytical result for 
the temperature dependent zero-bias conductance in Ref. 25 
and thus allows us to estimate (from the temperature at which 
I(V -> 0) = J /2) the Kondo temperature as T K /T « 0.24, 
in good agreement with the a priori estimate from the stan- 
dard formula 26 



T, 



K 



r \!/2 



-tvU/8T+tvT/2U 



(23) 



This formula is valid in the strong correlation regime and for 
U/T = 6 yields T K /T = 0.21. 

As V is increased the temperature dependence is weakened. 
At intermediate values of V, in the Coulomb blockade regime 
(V/T = 2), the current has little temperature dependence at 
low T. At large voltage bias (V/T = 4), correlation effects 
are already weakened due to the voltage, as is evident from the 
increase in I/Iq, and the almost perfect agreement with fourth 
order perturbation theory discussed in Section V. The current 
in this regime remains insensitive to temperature at low T. 



IV. RESULTS: VOLTAGE QUENCH 

A. Cutoff dependence 

An alternative procedure to calculate the steady state cur- 
rent of interacting quantum dots is to start from an interact- 
ing state in equilibrium (V = 0) and turn on the voltage at 
t = 0+. While this approach is computationally more ex- 
pensive and is restricted to nonzero temperatures, because it 
involves operators on the imaginary-time branch, it has the 
advantage that the Kondo resonance in the many-body den- 
sity of states is present already in the initial state. The Kondo 
resonance is built up during the evolution along the imaginary - 
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FIG. 6: Effect of the smoofhing-parameter v and cutoff lu c for 
f3F = 10. The top panel shows the interacting (U /F = 4) and non- 
interacting current for a quench to V/T — 0.25. The red line and red 
circles correspond to vT — 10 and cutoff uj c /T = 10, the blue lines 
and blue diamonds to vT = 3 and cutoff u) c /T = 10. A sharp band 
edge (large value of v) leads to oscillations in the current which make 
it difficult to estimate the steady state value. The bottom panel plots 
non-interacting and interacting currents for V/T — 0.5, vT = 3, 
f3T — 10 and indicated values of the cutoff. While the short-time 
behavior of the current is cutoff dependent, the steady state value is 
essentially cutoff-independent, as long as ui c > V. 



time branch, which does not add significantly to the sign prob- 
lem. One might expect that this V^-quench is particularly suit- 
able to study the small voltage regime, because turning on a 
small voltage will not change the spectral function dramati- 
cally. 

Since the voltage quench has not yet been discussed in the 
previous literature, we will now analyze the properties of the 
current after such a V-quench in some detail, and in particular 
its dependence on the band width (u> c ) and smoothness of the 
cutoff (v). The top panel of Fig. 6 shows the time evolution 
of the current in a model with U/T = (lines) and U/T = 
4 (symbols) if the voltage is suddenly increased to V/T = 
0.25 <C oj c /T = 10. In the initial state, the system is in 
equilibrium, with no current flowing through the dot. After the 
voltage bias is turned on, the current increases. In the model 
with hard band cutoff (^r = 10, red online) oscillations in the 
current appear which make it difficult to estimate the steady 
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FIG. 7: Voltage dependence of the current for U/T — 4 (top panel) 
and (U /T = 6) bottom panel. Solid lined show the non-interacting 
current and symbols the interacting current. As the voltage reaches 
V/T = 1 (U/T = 4) or V/T = 0.5 (U/T = 6), the interact- 
ing current overshoots, resulting in a slow convergence to the steady 
state. For smaller voltages, however, the steady state current can be 
computed on the basis of ^-quenches. All results are for /3T — 10, 
w o /r=10, and v/T = 3. 



state value. A smoother band cutoff (isT = 3, blue online) 
almost completely eliminates these oscillations. We will thus 
in the rest of this subsection show results for the "smoothing 
parameter" vT = 3. 

The lower panel of Fig. 6 illustrates the dependence of the 
current on the cutoff value oj c . While the short time behavior 
of the current depends strongly on the bandwidth, the steady 
state value shows little cutoff dependence as long as uj c is sub- 
stantially larger than the applied voltage bias. This is consis- 
tent with the observation for the U -quench in Ref. 11, where it 
was found that w c /T = 10 was enough to get accurate results 
up to V/T — 10. Hence, we will choose lu c /T = 10 for the 
rest of this paper. 



B. Voltage dependence 

In Fig. 7 we plot the non-interacting and interacting cur- 
rents for j3T = 10 and several values of the voltage bias. The 
top panel is for U/T = 4 and the bottom panel for U/T = 6. 




FIG. 8: Temperature effect on the current in the small voltage bias 
regime. Red circles, blue diamonds and black diamonds show the 
interacting current (U/T = 4) for V/T = 0.25 and V/T = 0.5, 
while the red, blue and black curves indicate the corresponding non- 
interacting currents. As the temperature is lowered, the interacting 
current increases towards the non-interacting value. 



In the small voltage regime (V/T < 0.5) the interacting cur- 
rent increases monotonically with time and eventually settles 
into a steady state within the accessible time window. The V- 
quench therefore allows us to measure accurate steady state 
currents for finite temperature in the small voltage regime. 
However, once the voltage becomes too big (see V/T = 1 
in the top panel of Fig. 7), the interacting current overshoots 
and only slowly settles into the steady state, making it im- 
possible to measure an accurate steady state value using this 
approach. However, as shown in the previous section, the 
simulation based on the U -quench provides accurate results 
once V/T > 0.5. The two simulation methods are therefore 
complementary in the sense that the F-quench works best for 
small voltage bias (V/T < 0.5) and the ?7-quench at larger 
voltage bias (V/T > 0.5). 



C. Temperature dependence 

The temperature dependence of the current after a V- 
quench in the low-voltage regime (V/T = 0.25 and 0.5) is 
shown in Fig. 8, which plots results for U/T = (lines) 
and U/T = 4 (symbols) for /3T = 5, 10 and 20. A rather 
strong temperature dependence is evident, in particular in the 
interacting current. This is consistent with the ?7-quench data 
shown in Fig. 5 and a consequence of the destruction of the 
Kondo resonance by temperature. Remarkably, a strong tem- 
perature dependence is observed even for T <C V, which 
means that the applied voltage does not effectively raise the 
temperature to a value of order V. In this voltage regime the 
non-zero voltage state therefore is not simply equivalent to a 
thermal state. The nature of the correlations which give rise 
to the temperature dependence are an interesting subject for 
further investigation. 
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FIG. 9: Current-voltage characteristics of the single-orbital Ander- 
son impurity model in the small voltage regime for U/T — 4, 6 at 
j3T = 10. The blue circles and black diamonds show V-quench re- 
sults. For comparison, we also plot finite temperature [/-quench data 
(U /r = 6, green crosses). 
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FIG. 10: Current- voltage characteristics of the single-orbital An- 
derson impurity model. The symbols show Monte Carlo data for 
U/T — 4, 6, 8 and 10, while the lines correspond to the fourth or- 
der perturbation calculation of Ref . 1 8 . The Monte Carlo results have 
been obtained by means of U -quenches at T — 0. Error bars are on 
the order of the symbol size. 



D. Comparison to the interaction quench 

In the F-quench calculations, for U/T = 6, we can access 
temperatures down to f3T w 20. At even lower temperatures, 
the perturbation order on the imaginary time branch becomes 
so large and the individual Monte Carlo updates so expensive 
that it is increasingly difficult to reach the very high statis- 
tical accuracy required for simulations with average signs of 
the order 10~ 3 . Since the problem of slow convergence in 
[/-quench calculations at small bias is considerably alleviated 
by finite temperature, it turns out that the accuracy of the lat- 
ter approach matches that of T^-quench calculations even at 
very small voltage bias (see U/T = 6 data in Fig. 9). For 
finite temperature simulations in the experimentally relevant 
temperature range {(3T > 20), the U -quench approach thus 
appears to be more powerful and sufficient to treat the entire 
voltage range. The good agreement between the [/-quench 
and V^-quench results in Fig. 9 furthermore shows that the 
steady state results obtained by the diagrammatic Monte Carlo 
method do not depend on the initial preparation of the system. 



V. I-V CHARACTERISTICS 

We now apply the machinery described in the previous sec- 
tion to compute the current-voltage characteristics of the An- 
derson impurity model at half-filling. The initial rise of the 
current at finite temperature (/3T = 10) is shown in Fig. 9. 
The blue circles and black diamonds have been obtained us- 
ing the V^-quench. The current-voltage characteristics in the 
V — ► limit becomes linear, although the slopes of the in- 
teracting and non-interacting models are not identical. This 
is the temperature effect on the Kondo resonance (particularly 
pronounced for large U) which was discussed in the previ- 
ous sections. As the temperature is lowered to zero, the ini- 
tial slope of the current approaches that of the non-interacting 



model. 

Figure 10 shows the T = result obtained using interac- 
tion quenches {lo c /T = vT = 10, essentially the wide band 
limit). The black curve shows the monotonic increase of the 
non-interacting current with increasing applied bias voltage. 
The red, blue and pink lines show the interacting current for 
U/T = 4, 6 and 8 predicted by fourth order perturbation 
theory. 18 Consistent with analytical arguments, 1415 the inter- 
acting current initially rises with the same slope as the non- 
interacting current, and reaches the non-interacting value also 
in the large-voltage limit. At intermediate values of V the ef- 
fect of interactions is to suppress the current (Coulomb block- 
ade). In fourth order perturbation theory, a hump appears in 
the I-V curve around V/T = 2 for U/T = 6 and 8. At even 
larger U (clearly outside the range of applicability) fourth or- 
der perturbation theory will presumably lead to a negative dif- 
ferential conductance at intermediate V. The appearance of 
this hump is related to the splitting of the Kondo resonance as 
discussed in Ref. 18. The Monte Carlo results for U/T = 4, 
6, 8, and 10 are shown by the red stars, blue circles, pink di- 
amonds, and orange triangles, respectively. Since these are 
U -quench results for T = 0, only V/T > 0.5 data are shown. 
In the large voltage regime (V/T > 4) the numerical results 
agree with the prediction from fourth order perturbation the- 
ory. Apparently, the fast decay of the Green functions for large 
voltage bias simplifies the diagram structure such that fourth 
order in £ is sufficient at V/T > 4. At intermediate voltages, 
1 IS V/T < 3, differences between the Monte Carlo data 
and fourth order perturbation theory appear. The essentially 
exact numerical data show no prominent hump feature near 
V/T = 2, and hence no negative differential conductance in 
the intermediate to strong correlation regime. The hump, and 
the associated splitting of the Kondo resonance, must there- 
fore be an artefact of fourth order perturbation theory. This is 
consistent with the conclusion reached in Ref. 1 1 on the basis 
of (less accurate) hybridization expansion results. The data in 
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Fig. 10 show that fourth order perturbation theory yields cor- 
rect results over the entire voltage range for U/T < 4. For 
larger interactions, and in particular around V/T « 2 more 
complicated self energy diagrams become important. 



of the U- or l/-quench. By slowly ramping up the interaction 
or voltage bias, it may be possible to avoid overshooting and 
thus observe a faster relaxation of the current into the steady 
state. 



VI. CONCLUSIONS 

We have discussed the implementation of the weak- 
coupling continuous-time Monte Carlo method on the L- 
shaped Keldysh contour and the application of this formal- 
ism to the study of transport through a quantum dot. Calcula- 
tions based on interaction quenches from the current carrying 
state at U = can be restricted to the real-time contours and 
provide accurate steady state currents for V/T > 0.5 and in- 
teraction strenghts U/T < 10, for arbitrary temperature and 
bandwidth. At finite temperature, convergence into the steady 
state is considerably faster, which allows access to the small 
voltage regime. As an alternative method, we have introduced 
calculations based on voltage quenches, which start from the 
interacting equilibrium state and which can be used to cal- 
culate the steady state current in the small voltage regime 
(V/T < 0.5) at finite temperature. Since the sign problem 
turns out to be essentially independent of the number of op- 
erators on the imaginary time branch, temperatures of order 
f3T = 10 can easily be dealt with. The V^-quench approach is 
however not more efficient than finite-temperature U -quench 
calculations. 

We have used the methods to accurately compute the 
current-voltage characteristics of the half-filled Anderson- 
impurity model in the intermediate-to-strong coupling regime. 
Comparison to fourth order perturbation theory showed that 
the latter fails at voltages around V/T « 2 for U/T > 4, 
but becomes accurate for V/T > 4. The splitting of the 
Kondo resonance predicted by low order perturbation theory 
is an artefact not present in the numerical data. The results 
presented in this paper show that diagrammatic Monte Carlo 
is one of the most powerful numerical tools for the study 
of non-equilibrium systems. The accuracy of the improved 
weak-coupling approach and its range of applicability rivals 
or surpasses other state-of-the-art numerical approaches such 
as time-dependent DMRG. 22 ' 23 For most practical purposes 
the numerical problem of calculating the steady-state current 
through a half-filled Anderson impurity model with symmet- 
rically applied voltage can be considered as solved. 

An interesting, and presumably straight-forward extension 
of our work would be the study of asymmetrically applied 
bias voltages. One of the optimizations of the Monte Carlo 
algorithm - the suppression of the odd perturbation orders 
- is however specific to the particle-hole symmetric model. 
Away from particle hole symmetry, odd perturbation orders 
contribute to the current and therefore must be considered in 
the simulation. This leads to an increase in the average per- 
turbation order and to a more severe sign problem, such that 
the accessible times and interaction strengths will be reduced. 
The optimal choice of the K x -parameters in the particle-hole 
asymmetric case is an open problem for future investigations. 
Another issue which should be considered is the optimal shape 
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APPENDIX A: NONINTERACTING GREEN FUNCTION 
FOR THE VOLTAGE QUENCH 

In this appendix we present the formalism needed for the 
voltage quench, evaluating noninteracting Green functions on 
the L-shaped contour with lead chemical potential /j, a equal to 
the equilibrium value /i(0) for times on the imaginary contour 
and with arbitrary time-dependence /i a (t) on the real time 
portions of the contour. The results of the paper correspond to 
H a (t)=n(0)±V/2. 

Because the voltage bias is time dependent, noninteracting 
Green functions cannot be expressed in the form of a Fourier 
transform, and instead they are computed numerically by the 
solution of their equations of motions in real (imaginary) time. 
A closed set of equations is obtained if one considers the non- 
interacting dot Green function [Eq. (13)], 

Go.a(t,0 = -i(Tcd^t)di(t')) , (Al) 
the hybridization of the dot to a single bath level 

G£„(i,f) = i(T e d a (t)a^(t')) , (A2) 
and the dot-decoupled Green function of a single bath state, 

= -i(7 C a° a (t) <i(0)o,v-=o. (A3) 

Here t, t' are arbitrary points on the real or imaginary por- 
tions of the contour, the time evolution is performed with 
U = but time-dependent voltage bias, and (-)o is the grand- 
canonical expectation value in the noninteracting initial state 
(at fj, = /x(0)). The contour-ordering operator Tq exchanges 
the product A(t)B(t') of two operators if t is earlier on the 
contour than t' (a minus sign is added if the exchange involves 
an odd number of Fermi operators). Equations of motions for 
the Green functions (Al) to (A3) are obtained from taking 
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time-derivatives and evaluation of the resulting commutators, 

[idf + e d )G , a (t,i?) =Y, G %A t > t 'K a ~ *c(t,0, ( A4 > 
p 

[id t , + e a pa - » a (t')} G a pa {t, U) = Go (V p a )*, (A5) 

[idv + e;„ - fi a (t')) 9p * a (t, t') = -6 c (t, t'). (A6) 

Note that when t = —it is on the imaginary branch of the 
contour, the time derivative is given by dt = id r . The con- 
tour delta function is given by 6c(t, t') = dtQc{t, t'), where 
Qc(t,t') = 1 if t is later on C than t' and zero otherwise. 
Eqs. (A4) to (A6) have a unique solution, provided that the 
contour Green functions satisfy an antiperiodic boundary con- 
dition on the contour C in both time arguments. 
Equation (A6) can be solved explicitly, 

g^(t,t')=g a (t,t';e; a ), (A7) 
g a (t,t';e) = -i[Q e (t,t') - f(e- fx(0))] 

x cxp J di[e~ ^ a (t)} — i J dt[e — fi a , (A8) 

where J* dt is along the contour. Furthermore, one can show 
from Eqs. (A5) and (A6) that the solution of Eq. (A5) is given 
by 

G^(M') = / d S G , a (t,s) (V p a )* g^(s,t'), (A9) 

where the integral runs over the whole contour. This expres- 
sion is inserted into Eq. (A4) in order to derive a single closed 
equation for Go, 

[id t , + e d ]Go,a(*, -Ida GoAt, »)M«> t') = -S c (t, t'). 

(A10) 

Here the sum over bath states has been condensed into the 



integral over the hybridization function (5) 

M*>0 =!>?(*. 0. (AH) 

OL 

A^(t,t')=J2K a \ 2 9^(t,t') (A12) 
de-T a (e)g a {t,t';e). (A13) 



In practice, we determine A from Eqs. (A13), (A8) and (5). 

Equation (A10) is an integrodifferential equation on the 
contour C. Its solution is equivalent to a boundary value prob- 
lem for the imaginary time component of the Green function 
and initial value problems for the components involving real 
time-arguments. The equation is solved numerically, using 
Langreth rules for the decoupling of real and imaginary time 
components [see Ref. 24]. 

The correlator (15) which enters Eq. (14) for the current is 
by definition given by 

A At,t') =iJ2 G P*(t,t')V p L v (A14) 
v 

Using Eqs. (A9) and (A12), this function can be obtained from 
the contour integral 



A . a (t,t') = / dsG .a{t,s)A^{ s ,t'). 



(A15) 



Note that the equations of motion (A5) and (A6) also 
hold in the interacting case, with Go replaced by G. Hence 
Eqs. (A14) and (A15) are still valid in the interacting case with 
the same replacement, and the interacting current can be ob- 
tained directly from the interacting dot Green function. This 
procedure is however equivalent to the approach which is used 
in the present paper, where the Green function is not measured 
and the current is obtained instead from Eq. (14). 
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